Different orderings in the narrow-band limit of the extended Hubbard model on the 
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We present the exact solution of a system of Fermi particles living on the sites of a Bethe lattice 
with coordination number z and interacting through on-site U and nearest-neighbor V interactions. 
This is a physical realization of the extended Hubbard model in the atomic limit. Within the 
Green's function and equations of motion formalism, we provide a comprehensive analysis of the 
model and we study the phase diagram at finite temperature in the whole model's parameter space, 
allowing for the on-site and nearest-neighbor interactions to be either repulsive or attractive. We 
find the existence of critical regions where charge ordering (V > 0) and phase separation (V < 0) 
are observed. This scenario is endorsed by the study of several thermodynamic quantities. 

PACS numbers: 71.10.Fd Lattice fermion models 71.10-w Theories and models of many-electron systems 
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I. INTRODUCTION 

In recent years, many theoretical as well as experi- 
mental investigations in condensed matter physics have 
been devoted to the study of low-dimensional strongly 
correlated electron systems where long-range Coulomb 
interactions play an important role. A crucial problem 
is to understand the effects of competing interactions 
and the corresponding phase transitions. One of the 
seminal models adopted to take into account long-range 
Coulomb interactions is the so-called extended Hubbard 
model (EHM), which, beside the on-site interaction U, 
contains also a nearest-neighbor interaction V . Although 
the EHM looks deceptively simple and in spite of a very 
intensive study, both analytical and numerical, there is 
no complete exact solution even in ID. These difficul- 
ties have led to the investigation of the so-called atomic 
limit of the extended Hubbard model (AL-EHM). Ac- 
cording to the conventional definition used in the liter- 
ature, the atomic limit stands for the classical limit of 
the model, where the hopping matrix elements are 
set to zero from the very beginning. The AL-EHM is 
exactly solvable in one dimension (see Ref. [1] and ref- 
erences therein), as well as when it is defined on the 
Bethe lattice [2, 3]. Here wc consider the Bethe lat- 
tice as the infinite version of a Cayley tree; thus we are 
not concerned with surface effects. The Bethe lattice, 
beside being an useful framework for studying the elec- 
tronic structure of amorphous and glassy solids, has a 
relevant role in condensed matter physics and statistical 
mechanics. Due to its peculiar structure, on the Bethe 
lattice it is possible to exactly solve several interesting 
physical problems involving interactions [4]. There are 
two special properties that make Bethe lattices particu- 
larly suited for theoretical investigations: the self-similar 
structure which may lead to recursive solutions and the 
absence of closed loops which restricts interference effects 



of quantum-mechanical particles in the case of nearest- 
neighbor coupling. Furthermore, Bethe and Bethe-like 
lattices have attracted a lot of interest because they usu- 
ally reflect essential features of systems even when con- 
ventional mean- field theories fail [5]. The reason is that 
such lattices are capable to take into account correla- 
tions which are usually lost in conventional mean-field 
calculations. Very recently, the Bethe lattice has been 
considered as the underlying lattice to study the suppres- 
sion of the paramagnetic metal- insulator transition in the 
Hubbard model at half-filling in the presence of nearest- 
neighbor and next-nearest-neighbor hoppings when the 
latter is increased with respect to the former [6, 7]. Fur- 
thermore, there exist in nature hyperbranched polymers 
with distinct regular molecular architecture (so-called 
dendrimers) which can be conveniently modelized by a 
Bethe-like Hamiltonian [8-10]. 

In this article, we study the AL-EHM on the Bethe 
lattice by means of the equations of motion approach 
[11]. For one-dimensional lattices, or more generally 
for lattices with no closed loops, classical fermionic and 
spin systems can be easily solved by means of the trans- 
fer matrix method [4]. However, this method is hardy 
implcmentable when more complex lattices are consid- 
ered. The Onsager solution for the two-dimensional Ising 
model is an emblematic example [12]. We feel that there 
is the necessity to foster alternative methods which can 
be used for a large class of lattices. In Ref. [13] we have 
shown that the equations of motion method provides an 
answer to this exigency. By using this method, local- 
ized fermionic systems and Ising-like models can be in 
principle solved for any underlying lattice. Indeed, one 
can find a set of eigenenergies and eigenoperators of the 
Hamiltonian which closes the hierarchy of the equations 
of motion. As a result, one can derive exact analytical 
expressions for the relevant Green's functions and corre- 
lation functions, which turn out to depend on a finite set 
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of parameters. The knowledge of these parameters is es- 
sential to obtain a solution of these models. In a series of 
articles [1-3, 13-17], we have analyzed several fermionic 
and spin systems and we have developed a self-consistent 
method which allows us to determine these parameters 
for the case of one-dimensional and Bcthc lattices. Some 
preliminary results for the AL-EHM on the Bethe lat- 
tice were presented in Refs. [2, 3], where we considered 
only the case of attractive intersite interactions (V < 0). 
Upon varying the temperature, we found a region of neg- 
ative compressibility, hinting at a transition from a ther- 
modynamically stable to an unstable phase, character- 
ized by phase separation. In this paper, we provide a 
comprehensive and systematic exact analysis of the AL- 
EHM on the Bethe lattice with coordination number z by 
considering relevant response and correlation functions 
as well as thermodynamic quantities. We consider both 
the cases V < and V > and cover a wide range of 
values of the parameters n, T and U (n is the particle 
density and T the temperature). The possibility for the 
parameters U and V to take positive as well as negative 
values - representing effective interaction couplings tak- 
ing into account also other interactions (for instance with 
phonons) - gives rise to a rich phenomenology and a vari- 
ety of phases. In particular, the ratio of the two compet- 
ing terms in the model Hamiltonian, i.e., the on-site and 
the intersite interactions, determines the different distri- 
butions of the electrons and the accessible phases, which 
we identified as charge ordered (CO) or phase separated 
(PS). Indeed, several studies of the EHM on regular lat- 
tices have pointed at the presence of states with phase 
separation [18-21] and of charge-ordered phases [20-25]. 
Furthermore, the introduction of an intersite interaction 
can mimic longer-ranged Coulomb interactions needed to 
describe effects observed in conducting polymers [26] . In 
this article, we report the phase diagram in the space 
(U, n, T) showing the existence of critical regions where 
charge ordering (V > 0) and phase separation (V < 0) 
are observed. These studies have not been reported in 
the literature and bring some more information on the 
properties of the Hubbard model. 

The plan of the paper is as follows. In Sec. II, we 
review and extend the analysis of the AL-EHM defined 
on a Bethe lattice with coordination number z leading to 
the computation of the Green's and correlation functions 
[2, 3] . In Sec. Ill, the attractive V case is reviewed in 
detail and the phase diagram is reported for a wide range 
of values of the external parameters n, T/ \V\, U/ \ V\. In 
Sec. IV, we investigate the case of repulsive intersite in- 
teractions. The phase diagram in the space n, T/V, U/V 
is derived and a long-range CO state is observed. The 
lattice is characterized by a inhomogeneous distribution 
of the particles in alternating shells. Relevant thermody- 
namic quantities - such as the specific heat, susceptibility, 
entropy - are also investigated as functions of the tem- 
perature, on-site potential and particle density. Finally, 
Sec. V is devoted to our concluding remarks while the 
appendix reports some relevant computational details. 



II. THE HAMILTONIAN AND THE 
EQUATIONS OF MOTION 

The theoretical framework leading to the exact solu- 
tion of the AL-EHM defined on the Bethe lattice has 
been already reported in Refs. [2, 3]. In this Section we 
shall review the analysis and extend it, considering also 
the breaking of translational invariance and the addition 
of an external magnetic field. 

In the extended Hubbard model, a nearest-neighbor 
interaction V is added to the original Hubbard Hamilto- 
nian, which contains only an on-site interaction U: 

<ij> 

U and V are the strengths of the local and intersite inter- 
actions, respectively; jj, is the chemical potential, n(i) = 
n-f(i) + ni(i) and D(i) = n-\{i)ni{i) = n (i) [n (i) — 1] /2 
are the particle density and double occupancy operators, 
respectively, at site i; ty is the hopping matrix. As usual, 
n„(i) = cl(i)c a (i) with a = {],[} where c a (i) (c+ (i)) is 
the fermionic annihilation (creation) operator of an elec- 
tron of spin a at site i, satisfying canonical anticommuta- 
tion relations. We use the Heisenberg picture: i = (i,t), 
where i stands for the lattice vector Ri. In the extreme 
narrow-band (atomic) limit the Hamiltonian (1) becomes 

H = -,*5>(i) + U + \ £ Vy n(i)n(j). (2) 

i i <ij> 

We shall study this model on a Bethe lattice with coor- 
dination number z. For this lattice, the Hamiltonian (2) 
can be conveniently rewritten as 

z 

H = - /m(0) + UD(0) + J2 # (p) • (3) 
P =i 

is the Hamiltonian of the p-th sub-tree rooted at 
the central site (0) and can be written as 

2-1 

ff(p) = -Lin( P ) + UD(p) + Vn(0)n(p)+Y^ H {p - m \ (4) 

m— 1 

Here (p) (p = 1, . . .z) are the nearest-neighbor sites of 
(0), also termed the first shell. i?( p ' m ) describes the m-th 
sub-tree rooted at the site (p); (p,m) (m = 1, ... z — 1) 
and (0) are the nearest-neighbors of the site (p). The 
process may be continued indefinitely. 

The equations of motion approach in the context of the 
composite operator method [11] - based on the choice of a 
convenient operatorial basis - provides us with the exact 
solution of the model. For our purposes, the suitable 
field operators are the Hubbard operators, £(i) = [1 — 
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n(i)]c(i) and r](i) = n(i)c(i), which satisfy the equations 
of motion: 



i^- t m = -K(*) + zVt(i)n a (i), 
i^ t nii) = {U- l i)r ] {i) + zVr l (i)n a [i). 



(5) 



In the following, for a generic operator $ (i) we shall use 
the notation <J> Q (?) = J2p=i *&{hP)/ z -> where (i,p) are 
the first nearest-neighbors of the site i. The Heisenberg 
equations (5) contain the higher-order nonlocal operators 
£(i)n a (i) and r](i)n a (i). By taking time derivatives of 
the latter, higher-order operators are generated. This 
process may be continued and an infinite hierarchy of 
field operators is created. However, since the number 
n(i) and the double occupancy D(i) operators satisfy the 
following algebra 



n p (i) — n(i) + a p D(i) 7 
DP (*) = !>(»), 
n" (i)D(i) = 2D(i) + a p D(i), 



(6) 



where p > 1 and a p = 2 P — 2, it is straightforward to 
establish the following recursion rule [13, 14]: 



(7) 



The coefficients 4m are rational numbers, satisfying the 



relations J2m=i Am = 1 and = <5 mjfe (k = 1, . . . , 2z) 
[1, 27]. The recursion relation (7) allows one to close the 
hierarchy of equations of motion. 



A. Eigenoperators and eigenvalues 

One may define the composite field operator 



(8) 



where 



( \ 

V(i)[n a (i)} 

(9) 

By using the recursion rule (7), one can show that the 
fields and ^ v \i) are eigenoperators of the Hamil- 

tonian (3) [2, 3]: 



£(i)[n"(i)] 



(10) 



where and are the (2z + 1) x (2z + 1) energy 
matrices 
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[/ - A* + zVAg z+1) J 



whose eigenvalues, E„ and Em , are given by: with to = 1,...,2z + 1. The Hamiltonian has now been 

formally solved since, for any coordination number of the 
E$ = —fx + (to — l)V, underlying Bethe lattice, one has found a closed set of 

EW = -» + U+(m-l)V, U,>) 
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cigcnoperators and eigenenergies. As a result, one may 
compute observable quantities. This will be done in the 
next section by using the formalism of Green's functions 
(GF). 



B. Retarded Green's functions and correlation 
functions 

The knowledge of a set of eigenoperators and eigenen- 
ergies of the Hamiltonian allows one to find an exact 
expression for the retarded Green's function 

G«(t - *') = 9{t - 0<{V> (s) (M)^ (s)t (M')})> (14) 

and, consequently, for the correlation function 

CW(t-t') = (V (s) (i,^ (s)t (i,i'))- (15) 

In the above equations, s = £,77 and (• • • ) denotes the 
quantum-statistical average over the grand canonical en- 
semble. It is not difficult to show that, for fermionic 
operators, only the on-site correlations are non-zero. It 
is not difficult to show that the retarded GF and satisfies 
the equation 



lj — e 



(16) 



where G^(lu) is the Fourier transform of G^ s \t — t') 

and /( s ) = ({V> (s) 00,-(/> (s)t (i)}) is the (2z + 1) x (2z + 1) 
normalization matrix. The solution of Eq. (16) is [11]: 

2z + 1 a (s,m) 



(17) 



uj - E$ + iS 

Similarly, the correlation function satisfies the equation 



1 + tanh — 



Im 



(18) 



where C^(uj) is the Fourier transform of C^(t — t'), 
whose solution is 



2z + l 



(s) 

I a h are instead the elements of the normalization matrix 



l( s \ Calculations show that 
matrix f2 given by 



fl(i) = 0, with the 



p,k 



1(A) 



1 



2z+l-p 



k = 1, p = 1, 
k = l, p^l, (21) 
fc^l. 



The matrix elements of the normalization matrices 1^ 
have the expressions 

j(0 — K (™+ m - 2 ) _ ^(n+m-2) j(rj) _ ^(n+m-2) ^2) 

where the correlators and A^ p ^ are defined as 

« W = (KW] P >, AW = l(n(i)[n«(i)F>. (23) 

By exploiting the recursion relation (7), it is not difficult 
to show that also and obey similar recursion 
relations 



2z 



2z 



? (P) 



E ^ ) « (m) . aw = e^ } a ( 



m) 



(24) 



m— 1 



m— 1 



limiting their computation to the first 2z correlators [2]. 
At this stage, the knowledge of the GFs and of the CFs 
is not yet achieved since they depend on the {<r( s ' TO )} 
which, in turn, depend on the normalization matrix el- 
ements: there are 2z parameters to determine. To find 
these parameters, we shall exploit the Pauli principle and 
impose pertinent boundary conditions for obtaining a set 
of self-consistent equations. 

One first chooses an arbitrary site, say i, then splits the 
Hamiltonian (3) in the sum of two terms: H = + 
Hj l \ where 



H, 



P = -n[n(i) + zn a {i)] + U[D(i) + zD a (i)] 

p—l m—l 

= zVn(i)n a (i), 



(25) 



C (s) (lu) = 7T E °" (5 ' m) T ™ S (u - E m)- (I 9 ) and introduce the if^ -representation: the statistical av- 



m— 1 



erage of any operator O can be expressed as 



In the above equations, T$ = 1 + tanh (ftE^ /2), f3 = 

l/k B T and the E^ are given in Eq. (13). The spectral 

density matrices o^ n ^ can be computed by means of the 
formula [11]: 



(O) 



(e-^)o,i 



(26) 



T (s>") 



2z+l 



V \rt s 



As) 

1 cb ■ 



(20) 



c=l 



In Eq. (20), fi (s) is the (2z + 1) x (2z + 1) matrix whose 
columns are the eigenvectors of the energy matrix e^- s \ 



The symbol (• • • ) ,i stands for the thermal average with 

respect to the reduced Hamiltonian H^: i.e., (• • • ) ,i = 

Tr{- ■ ■ e-l 3H o ) }/Tr{e-t !H ° ) }. Equation (26) allows us to 
express the thermal averages with respect to the complete 
Hamiltonian H in terms of thermal averages with respect 
to the reduced Hamiltonian Hq, which describes a system 
where the original lattice has been reduced to the site i 
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and to z unconnected sublattices. As a consequence, in 
the Tip-representation, correlation functions connecting 
sites belonging to disconnected sublattices can be decou- 
pled. Let us consider the correlation functions 



C[ s i = (s{i)s^{i) [n a (it- 1 ), 



(27) 



where s — and k = 1, ...,2z + 1. By means of Eq. 
(26), they can be written as: 



(J) _ Mi) 8 t(j) [n°(i)] fc - 1 e-^)o, 1 



(e-^)o,i 



(28) 



The Pauli principle leads to the following algebraic rela- 
tions 



(, (i)n(i) = 0, Ty(i)n{i) = ry{i), 

e t (*)i>(*) = o, rf(i)D(i) = 0, 



(29) 



from which one has ^(i)e ^ Hl = ^(i), and 
rf{i)e-P H ' = rf{i)e- z P Vna ^. In' the H - 
representation, the correlation functions can be rewritten 
as: 



r (o _ (e(»)e t (»))o, i (K(»)] fc - 1 )o4 

hk (e-^) ,i 

= fa(0>? t (»))o,i([n a (»)] fc - 1 e-^ Vn ° (i) )o,i 
(e-^)o,i 



(30) 



In the iio-represcntation, the Hubbard operators obey 
to simple equations of motion: [£(i),ifo] = — an( i 
[77(1), i?o] = —(a* — ^) ^W- Thus, it is easy to show that 
the equal time CF's can be expressed as: 



1 



1 + 2e^ + eJ^-U) 



l-Bi + Ba, 



+ 2e ^ + e«2^c/) = 2 (Bl _ 2B2) ' 

(31) 



where: 



52 - W*)>«m - TTWT- g - 



(32) 



In Eq. (32), we have defined / = e^, g = e^- u \ and 
we have used the identities 



(33) 



Upon inserting Eqs. (31) into Eqs. (30) and by taking 
k = 1, one finds: 

(£) _ l-B 1+ B 2 

= (g 1 -2g 2 )(e-^"°W) , 1 W 
2(e-^)o,i 



It is not difficult to show that the averages in the above 
equations can be expressed as: 

{e- 0H ^ )o,i = 1 + Bi (f? - 1) + B 2 (1 - 2J? + G't) , 

(e -z/JVn-(0) )0)i = 

(35) 

where 

Fi = 1 + aXi + a 2 Yi, G t = 1 + dX, + d 2 Yi. (36) 
In the above equations we have defined a = K — 1, 



K 1 - 1, with if 



-pv. 



> *p (P — 1; • • • 7 Z ) i S an 



arbitrary neighboring site of i. X, and Yi are two param- 
eters defined as: 



A 4 = (" Q (*)}o,i = -^>(* P ))o,i 



P =i 



^ = (^W) ,i = i^(^ P ))o,i. 



(37) 



P =i 



Xi and Yj are parameters of seminal importance since 
all correlators and fundamental properties of the system 
under study can be expressed in terms of them. Relevant 
physical quantities, such as the mean value of the particle 
density and doubly occupancy, and the charge correlators 
and A( fe ) (23) can be easily computed. After lengthy 
but straightforward calculations, one finds: 



(n(i)) = 
(D(i)) = 



2fF* + 2gGl 
l + 2fF?+gGl 

gGj 

\ + 2fF?+gG\ 



(38) 



and 



(n(i p )> 



[2fK {X i + 2aY i )F t 



z-l 



(D(i P )) = 



\ + 2fF?+gG' i 
+ gK 2 {X l + 2dY l )Gr 1 +X t ] , 
Y l (l + 2fK 2 Fr 1 +gK i Gr 1 ) 



l + 2fF*+gGf 



(39) 



The parameters Xi and Yi will be fixed by using bound- 
ary conditions, which will be different according to the 
sign of the intersite potential V. Therefore, we shall con- 
sider separately the cases V < and V > 0. For both 
cases we shall study, in the next sections, relevant ther- 
modynamic quantities and response functions, such as 
the internal energy, the specific heat, the charge and spin 
susceptibilities and the entropy. The internal energy E 
can be computed as the thermal average of the Hamil- 
tonian (2) and it is given by E = UD + zV\ w . The 
specific heat is then directly given by C = dE/dT. The 
charge susceptibility %c can be computed by means of 
thermodynamics through the formula 



Xc = Nn 2 



1 dn 



(40) 



6 
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FIG. 1: Distribution of the particles in the Bethe lattice at 
n = 0.5 and T = 0: (a) V < and U > U PS ; (6) V < 
and U ;$ t/p,g. White, grey and black circles denote empty, 
arbitrary spin singly occupied and double occupied sites, re- 
spectively. 



In the above equation, N is the number of sites and 
n = '^2 i (n(i)} /N is the particle number per site. The spin 
magnetic susceptibility \s can be computed by introduc- 
ing an external magnetic field h, taking the derivative of 
the magnetization m = (n^(i) — n^(i)) with respect to h 
and letting h going to zero: 



Xs 



dm 
~dh 



(41) 



h=0 



The addition of a homogeneous magnetic field does not 
dramatically modify the framework of calculation given 
in this section, once one has taken into account the break- 
down of the spin rotational invariance. Some details of 
the calculations are given in the appendix. 



mine and Yi as functions of the chemical potential 

X = 2/(1 - X - dY)F z - 1 + g[2 + (d- 1)X - 2dY]G z ~\ 

(42a) 

Y = g[l + dX - (2d + l^G*" 1 - 2fK 2 YF z -\ (42b) 

In the above equations, we dropped the index i be- 
cause all the sites are equivalent. Since experimentally, 
by varying the doping, it is possible to tune the density 
in a controlled way, here we shall fix the particle density 
n = (n(i)); the chemical potential will be determined by 
the system itself, according to the values of the external 
parameters, by means of the following equation: 



(X - 2Y) F + 2YG 



1 — X + Y + (X — 2Y) F + YG 



(43) 



Equations (42) and (43) constitute a system of coupled 
equations allowing one to ascertain the three parameters 
/i, X and Y in terms of the external parameters of the 
model n, U, V, T, and z. Once these quantities are 
known, all the properties of the model can be computed. 
As an example, the double occupancy D = (D(i)) and 
the nearest-neighbor charge correlation function A*- 1 ) = 
(n(i)n a («)]) / '2, in terms of X and Y, are given by: 



D 



YG 



A (1 



1— X + Y + (X- 2Y) F + YG ' 
K (X + 2aY) (X - 2Y) + 2K 2 Y (X 



2dY) 



2 (1 - X + Y) + 2 (X - 2Y) F + 2YG 



(44) 



III. ATTRACTIVE INTERSITE POTENTIAL 

In this Section, we review the case of attractive inter- 
site potential. For V < 0, the AL-EHM on the Bethe 
lattice exhibits a phase separation at low temperatures 
[3]. This phenomenon is characterized by a macroscopi- 
cally inhomogeneous ground state where different spatial 
regions have different average particle densities. As it has 
been evidenced in Ref . [3] , there exists a critical tempera- 
ture T c below which the system loses translational invari- 
ance and phase separation occurs. Moreover, for T < T c 
two types of phase separated configurations are possible 
according to the strength of the on-site potential U : clus- 
ters of singly occupied sites or clusters of doubly occupied 
sites. This is qualitatively illustrated in Fig. 1, where we 
report just one possible configuration for n = 0.5 and 
T = 0. There is a critical value of the on-site poten- 
tial Ups, depending on the coordination number z and 
on filling n, separating the two types of phase separated 
configurations. In the high temperature regime, one may 
safely assume that the system is in the translational in- 
variant phase. Thus, one requires: (n(i)) = {n a (i)) and 
(D(i)) = (D a (i)), Mi. As a consequence, from Eqs. (38) 
and (39), we obtain two equations allowing us to deter- 



A. The solution at half filling 

As it has been evidenced in Ref. [2], of particular 
interest is the case of half filling due to the possibility 
to analytically reveal the spontaneous breakdown of the 
particle-hole symmetry. Because of this invariance prop- 
erty of the Hamiltonian (3), at n = 1 the chemical poten- 
tial does not depend on the temperature and takes the 
constant value /i = U/2 + zV. Upon substituting this 
value of /i in Eqs. (42), one finds that these equations 
admit the following solution 

X = l-dY, (45a) 

with Y determined by the equation 

Y(l + K 2 ) + 2e f}U ' 2 KY (l - a 2 Y) z ^ -1 = 0. (45b) 

It can be shown that this solution gives 

n = l, 

° = 2 + 2eW 2 (l - a?YY ' (46) 
A« = ± + dYD. 



u 



(a) 



FIG. 2: The transition temperature T c as a function of the 
on-site potential [/ for V = — 1, n = 1 and different values of 
the coordination number z. 



That is, solution (45) satisfies the particle-hole symmetry. 
On the other hand, if one perturbs the solution (45a), by 
setting for example X = 1 — dY + 5, it is straightforward 
to show that a particle-hole symmetry breaking solution 
exists for temperatures lower than a critical temperature 



T c , determined by: 



U/\V\+2 

2K^^[z + K(z-2)f 



+ (K + l) z -\z - ly-^z - K 2 (z - 2)] = 0. 



(47) 



Numerical calculations show that Eq. (47) admits a 
solution only for attractive intersite interactions. It is 
interesting to consider Eq. (47) in two extremal limits, 
namely U j \V\ — > -co and T c — > 0. In the former limit, 
corresponding to the case of vanishing intersite potential, 
one finds the same critical temperature exhibited by a 
system of spinless fermions living on the sites of a Bethe 
lattice [2]: 



In 



z-2 



(48) 



Since the spinless fermion model can be mapped into the 
spin- 1/2 Ising model, our result for the critical tempera- 
ture agrees with the one previously found in the literature 
[4, 15]. In the second limit, Eq. (47) becomes 



(U/2-\V\) = k B T c \n 



(*-!)'• 
2(z - 2) 2 



(49) 



That is, U = 2\V\ is the critical value of the on-site 
potential at which a quantum phase transition occurs. 
The results obtained from Eq. (47) are displayed in Fig. 
2. One observes that, at fixed coordination number, by 
increasing U from large negative values, the critical tem- 
perature decreases, and the lower the coordination num- 
ber, the lower the critical temperature. An interesting 
feature of the phase diagram is that for U > 2\V\, the 
critical temperature exhibits a reentrant behavior. The 
width of the reentrance increases with z, and the turning 
point is at a critical value U c , which depends linearly on 




(b) 



FIG. 3: The charge susceptibility \c for V = — 1 and z = 3 as 
a function of the temperature T for several values of U and 
(a) n = 0.5; (b) n = 1. 



the coordination number via the law: J7 C /|V| = a + bz 
(for z > 3), where a w 0.69 and b w 0.54. For z = 3 one 
has U c « 2.3|V|. 



B. Phase diagram and local properties 

In this Subsection, we review the case of arbitrary fill- 
ing analyzed in Ref. [3]. By solving the set of equations 
(42) and (43), it is possible to derive the phase diagram 
and various local properties in terms of the external pa- 
rameters n, U, and T, taking \V\ as the unit of energy. 
An important quantity useful for studying the critical be- 
havior of the system is the susceptibility. In Figs. 3, we 
plot the charge susceptibility as a function of the temper- 
ature for n — 0.5 and n = 1, and for several values of the 
on-site potential. From Fig. 3a, one observes that there is 
a critical temperature - depending on n and U - at which 
the charge susceptibility \c diverges, for both attractive 
and repulsive U. Below T c , \ c becomes negative (not 
reported in the graphs). As a result, one immediately 
infers that there exists a region where the system is ther- 
modynamically unstable. At fixed U, the instability is 
observed in a particle density region An = m < n < n 2 , 
whose width varies with the temperature, and vanishes 
for T > T c [3] . For values of the particle density less than 
half filling, the divergence of the charge susceptibility is 
observed for all values of U. As it is shown in Fig. 3b, for 
n = 1 and for U > U c , no phase transition is observed: 
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FIG. 4: The spin magnetic susceptibility \ s as a function of 
the temperature for V = — 1 and 2 = 3, and for (a) (7 = — 1 
and 7i = 0.25, 0.5, 0.75, 1; (b)U = 3 and n = l. 



Xc is well defined for all values of T and vanishes in the 
limit T — > 0. It can be shown that, in the opposite limit 
T — > oo, Xc tends to a constant value which does not 
depend on U but only on n according to the law 

f n 
lim Y r = n 1 

In Figs. 4, we plot the spin susceptibility Xs as a function 
of the temperature. For attractive intersite interaction, 
at T = T c one does not observe substantial changes in the 
behavior of Xs- For T > T c , the system is in the high- 
temperature regime and Xs follows a Curie law. As it 
is evident from Fig. 5b, the only low-temperature region 
accessible to our investigation is when n = 1 and U > U c . 
In Fig. 4b, we show the spin magnetic susceptibility Xs 
as a function of T for n = 1 and U = 3. For U > U c , Xs 
is rather insensitive to the value of U, and a large peak 
is observed at low temperatures. The behavior of the 
spin susceptibility for repulsive intersite interactions is 
more relevant, as we will show in the next section, since 
it signals the occurrence of a CO phase. 

In Fig. 5a, we report the phase diagram in the 3D 
space (U,n,T). The critical temperature T c and the 
width of the instability region An increase by decreasing 
U: an attractive on-site potential will favor phase sepa- 
ration and in particular the clustering of doubly occupied 
sites. As a function of the particle density, the critical 
temperature shows a lobe-like behavior: it increases by 
increasing n up to half filling, where it has a maximum; 
further augmenting n, it decreases vanishing at n — 2. A 




FIG. 5: (a) (Color online) The critical temperature T c as a 
function of the particle density n and of the on-site potential 
U for V = — 1 and z — 3. (b) The critical temperature T c as 
a function of the on-site potential U for V = — 1, z = 3 and 
for n = 0.25, 0.5, 0.75, 1. 



different behavior is observed when U > 2 \ V\: the lobe 
splits in two and T c increases with n up to quarter fill- 
ing, then decreases with a minimum at half filling. The 
critical temperature at n = 1 is finite only for U < U c . 
From Fig. 5b, one immediately notices that, for larger 
values of the on-site potential, there is no transition at 
half filling, as it was noticed in Ref. [3]. 

As it has been already pointed out, below the critical 
temperature T c the system is thermodynamically unsta- 
ble. However, the behavior of relevant thermodynamic 
quantities below T c unveils the existence of a critical 
value of the on-site potential separating the two types 
of phase separated configurations sketched in Fig. 1 [3]. 
Interestingly, this critical U has the same value of Ups, 
turning point in the T c -curve at half filling (see Fig. 2). 
In Figs. 6, we plot the double occupancy D, the short- 
range correlation function A^ 1 ^ and the internal energy E 
as functions of U for n = 0.75 and different values of the 
temperature. One observes that D and A^ 1 -* exhibit two 
plateaus at low temperatures. For z = 3, in the limit 
T — > 0, there is a discontinuity around Ups ~ 2.3|V|. 
For U > Ups, the double occupancy vanishes, whereas 
AW tends to n/2. The repulsion between the electrons 
on the same site and the concomitant nearest-neighbor 
attraction, leads to a scenario where the electrons tend 
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FIG. 6: (a) The short-range correlation function A' 1 ', (b) 
the double occupancy D and (c) the internal energy E as 
functions of U for V = -1, z = 3, n = 0.75 and for T = 0.2, 
0.3, 0.4, 0.5. 



behavior of D, A^ and C as functions of U. 

IV. REPULSIVE INTERSITE POTENTIAL 



to cluster together occupying neighboring sites. At half 
filling all sites are singly occupied, whereas for n < 1 
one observes two separated clusters of filled and empty 
sites. When U < Ups, one observes a dramatic increase 
(step-like as T — > 0) of the double occupancy and of 
the short-range correlation function, namely: D — > n/2 
and A^ 1 ) — ► n. As a consequence, the sites are dou- 
bly occupied, and A' 1 ) = n indicates that the doublons 
(the charge carriers of doubly occupied sites) tend to oc- 
cupy nearest-neighbor sites arranging to form large do- 
mains occupied, leaving the rest of the lattice empty. For 
U < Ups one also observes a dramatic decrease of the 
internal energy E, with a discontinuity as T — ► T c around 
U ~ Ups- In Figs. 7, we plot the specific heat C as a 
function of the temperature T at n = 0.75 and for several 
values of the on-site potential U. For attractive U, the 
specific heat presents one peak whose height decreases 
by decreasing U (see Fig. 7a). In Fig. 7b, one observes 
that for repulsive U the height of the peak increases by 
increasing U moving to lower temperatures. In the limit 
U — > Ups, the peak becomes sharper and a divergence 
is observed at U = Upg, confirming the emergence of 
a critical point. The same analysis carried out for dif- 
ferent values of the particle density, evidences a similar 



A repulsive intersite interaction disfavors the occupa- 
tion of neighboring sites. At low temperatures, this may 
lead to a CO phase characterized by a distribution of the 
electrons in alternating shells. In order to capture this 
phase, we shall solve the self consistent equations (38) 
and (39) by releasing the translational invariancc require- 
ment. To this end, one can divide the lattice into two 
sublattices: A contains the central point (0) and the even 
shells, the sublattice B contains the odd shells. Then, one 
requires the following boundary condition (BC) to hold: 

i 

Let us take two distinct sites i E A and j e B. We 
require that the expectation values of the particle den- 
sity and of the double occupancy operators at the site i 
are equal to the ones of the neighboring sites of j and 
viceversa, namely: 

(n(i)) = (n a V)), (n(j)) = {n a (i)), 
(D(i)) = (D"(j)), (D(j)) = (D"(i)). 
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Thus, by means of Eqs. 
equations 

2(fF A +gG A ) 



(38) and (39), one finds two 



1 



l + 2fFX+gG* A l + 2fF B+g G\ 



2fK(X B + 2aY B )F z B 1 + gK 2 (X B + 2dY B )G z B 1 ] , 
gG\ _ Y B (1 + 2fK 2 F z B - 1 + gK i G z B 1 ) 



l + 2fF%+gG* 



1 + 2fF* + gG z 



(52) 



plus two more equations which can be obtained by sub- 
stituting A <-» B. The coefhcients subscripts pertain to 
sites belonging to the two different sublattices. In order 
to close the set of self-consistent equations one needs one 
more equation. By using the BC (50), one can fix the 
particle density as 



fF% + gG* 



fF B + gG\ 



\ + 2fF% + gG A ! + 2fF z B +gG\ 



(53) 



As a result, one can now determine all the unknown pa- 
rameters X A , X B , Y A , Y B , and /i. By means of the same 
analysis employed in the previous section, one can ex- 
press the local correlators and A^ p ^ in terms of these 
parameters, and eventually compute all the properties of 
the system. In particular, one finds: 



D 



1 



(D A 



^D b ) 



2(1 + 2/F 
fK(X A + 



9+9G z A 
2aY A )F z A - 1 



2(1 + 2fF z B +gG z B Y 
gK 2 (X A + 2dY A )G z 



l + 2fF z B +gG\ 



(54) 



By varying the external parameters U, n and T one has 
the tools to completely characterize the phase diagram 
pertinent to the case of repulsive intersite interactions. 
In the following we shall set V = 1. 



A. Phase diagram 

In this section, we derive the phase diagram: by nu- 
merically solving the set of equations (52) and (53) we 
find a region of the (U,n,T) 3D space characterized by 
a spontaneous breakdown of the translational invariance. 
In this region, shown in Fig. 8a for the case z = 3, the 
population of the two sublattices A and B is not equiv- 
alent: the system has entered a finite temperature long- 
range CO phase. Upon decreasing the temperature, the 
distribution of the electrons becomes more and more in- 
homogeneous. To better understand the typology of the 
critical region we take sections of the 3D structure and 
study the 2D phase diagrams at constant n in Fig. 8b 
and at constant U in Fig. 9. 

In the plane (U, T), the critical temperature line shows 
different behaviors according to the value of the particle 
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FIG. 8: (Color online) (a) Phase diagram in the space 
(U, n, T) for V = 1 and z = 3. (b) Phase diagram in the 
plane (U, T) for V = 1 and z = 3 and several values of n. 



density. As illustrated in Fig. 8b for the case z = 3, one 
encounters the following situation, (i) For n < 1/z there 
is no CO phase for all values of U; (ii) For 1/z < n < n p 
a CO phase is observed only for U > (n p = 0.46 for 
z = 3); the critical temperature T c increases by increas- 
ing U and tends to a constant (depending on the value 
of n) in the limit U — > oo. (iii) For n p < n < 2/z there is 
a CO phase for both attractive and repulsive U, with a 
reentrant behavior for U < 0. (iv) For n > 2/z there is no 
reentrant behavior: from a constant value at large neg- 
ative U, the critical temperature decreases by increasing 
U and vanishes at a certain value of U . An interesting 
feature is the presence of a crossing point in the criti- 
cal temperature curves around U — 2.3V for z = 3. In 
the plane (n,T), at fixed U, the CO phase is observed 
in the interval m < n < ni\ the width A = — n\ 
varies with the temperature, following different laws ac- 
cording to the value of U. At T = a complete CO state 
is established in the regions 2/z < n < 2(z — l)/z for 
U < and 1/z < n < (2z - l)/z for U > 0. For U < 
(see Fig. 9a), An first increases with T, then decreases 
vanishing at n = 1 , where the maximum critical temper- 
ature is reached; a reentrant behavior characterizes this 
region. As it is shown in Fig. 9b, the value U = is 
a singular point; at T = the CO phase exists in the 
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FIG. 9: Phase diagram in the plane (T, n) for V = 1: Figs. 
9a, 9c, 9d 2 = 3 and several values of C/ ; Fig. 9b, U = and 
2 = 3,4,5. 
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interval n p < n < 2 — n p . For < U < 2 (see Fig. 9c), 
An decreases with T and no reentrant phase is observed; 
as in the attractive case, the phase diagram presents a 
single lobe structure centered at n = 1. For U > 2 (see 
Fig. 9d), one observes the formation of two lobes cen- 
tered around n = 0.5 and n — 1.5 and the corresponding 
decreasing of the central lobe centered at n = 1, which 
disappears at U = (U~o ~ 3.7 for z = 3). For U > 
the two lobes are separated; at T = the CO phase is 
observed in two regions centered at n = 0.5 and n = 1.5. 
It is worthwhile to notice that for U < 2V and for all 
values of T , An increases by increasing the coordination 
number z; in the limit z — > oo, An — > 2 and the CO phase 
is observed for all values of n. Besides the translational- 
invariance broken solution in the 3D region illustrated in 
Fig. 8a, the set of equations (52) and (53) admits also a 
homogeneous solution. In order to determine which so- 
lution is energetically favored one has to look at the free 
energy. The Helmholtz free energy F can be computed 
by means of the formula 



n 

F(T, n) = j n(T,n')dn'. 
o 



Upon defining the difference AF = Fh om — Fco - where 
Fhom is the free energy of the homogeneous phase and 
Fco the one of the CO phase, respectively - one finds 
that below the transition temperature the CO phase is 
energetically favored since AF > 0. In particular, AF 
smoothly vanishes for T — > T c , signalling a second-order 
phase transition. 



FIG. 10: (a) The chemical potential )i as a function of the 
temperature for V = 1, z = 3, n = 0.75 and different on-site 
interactions, (b) The chemical potential fi as a function of 
the particle density n for z = 3, T = 0.2 and various on-site 
interactions. 



B. Thermodynamic properties 

In this section, we shall determine several thermody- 
namic quantities whose behaviors support the scenario 
depicted in the previous section. The behavior of the 
chemical potential as a function of the temperature and 
of the particle density n is reported in Figs. 10a and 10b, 
respectively. One can immediately notice that, as a func- 
tion of the temperature, [i presents a cuspid at T = T c : 
for T > T c (T < T c ) /i is a decreasing (increasing) func- 
tion of T. When plotted at fixed temperature, the chemi- 
cal potential is always an increasing function of n, hinting 
at a thermodynamically stable system, for both the CO 
and homogeneous phases. In the limit T — > 0, /i shows 
two plateaus for U < 2V in the range < n < 2; when 
U > 2V each plateau splits in two sub-plateaus. 

A full comprehension of the phase diagram and of the 
distribution of the particles on the sites of the Bcthe lat- 
tice can be achieved by a detailed investigation of the par- 
ticle density and of the double occupancy. In particular, 
the study of the different contributions coming from the 
two sublattices, reveals the onset of CO states, the reen- 
trant behavior, and also the distribution of the electrons 
in the shells of the Bethe lattice. In Figs. 11a and lib, 
we plot the sublattices variables ua and ng as functions 
of the temperature for n = 0.5, n = 0.75, and several val- 
ues of U. At high temperatures the particles distribute 
homogeneously in the entire lattice: tia — ns — n. Upon 
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FIG. 11: The sublattices variables ua and n B as functions of 
the temperature for V = 1, z = 3, several values of {7 and: 
(a) n = 0.5, (b) n = 0.75. 



FIG. 12: The sublattices variables Da and Db as functions 
of the temperature as functions of the temperature for V = 1, 
z = 3, several values of £7 and: (a) n = 0.5, (b) n = 0.75. 



lowering the temperature, one finds that at the critical 
temperature T c — T c (n, U) a CO phase is established; the 
particles tend to fill one sublattice (B) and to empty the 
other sublattice (A). By further lowering the tempera- 
ture, the behavior is different according to the values of n 
and U. For n = 0.5, one observes the following situation: 
(i) for U negative and close to zero the phase diagram 
exhibits a reentrant behavior (see Fig. 8b); correspond- 
ingly, as shown in Fig. 11a, there is another critical tem- 
perature at which the translational invariance is restored 
and the two sublattices are again equally populated; (ii) 
for {7 = 0, the anisotropy of the filling increases and, at 
T = 0, one finds ua = 2n/5 and lis = 8n/5 for z = 3; 
(iii) for U > a total charge order is established in the 
limit T — > 0: all the particles reside in one sublattice (B) 
whilst the other is completely empty. For n = 0.75 and 
n = 1 (not shown), there is no reentrant phase; by low- 
ering T, ua («.b) decreases (increases) and tends to zero 
(2n) at T = (for n = 0.75 and U < charge ordering is 
not complete, since a small fraction of particles is found 
in the sublattice A, also at T = 0). 

In Figs. 12a and 12b, we plot the sublattice double 
occupancies Da and Db as functions of the tempera- 
ture for n = 0.5 and n = 0.75 and several values of U. 
At high temperatures the system is in a homogeneous 
phase and, correspondingly Da = Db = D. Upon low- 
ering the temperature, one finds that at T c there is a 
phase transition to a CO state: Db increases while Da 
decreases. By further lowering the temperature, one ob- 



serves different behaviors according to the values of n and 
U. A reentrant phase characterized by a second critical 
temperature at which the system returns to the homoge- 
neous phase where the double occupancy tends to n/2 in 
the two sublattices (observed for n — 0.5 and U negative 
and small). The CO persists in the limit T — > with: 
(i) Da and Db both vanishing at T = (observed for 
n = 0.5 and U > 0); (ii) Da (Db) decreases (increases) 
and tends to a finite value (observed for n = 0.75 and 
U < 0); (iii) Da (Db) decreases (increases) and tends to 
zero (n) (observed for n = 0.75, U > and n = 1, VU). 

Since the specific heat exhibits a very rich structure in 
correspondence to the critical lines in the phase diagram 
shown in Fig. 8a, we shall investigate its behavior for 
different values of the particle density n. The possible 
excitations of the ground state are creation and annihila- 
tion of singly occupied or doubly occupied states, induced 
by the Hubbard operators ip^ and tj}^\ respectively [1]. 
The corresponding transition energies are given in Eqs. 
(13) and the high and low temperature peaks exhibited 
by the specific heat are due to these transitions. One 
may distinguish between them by looking at the posi- 
tion of the peaks, i.e., if the position changes or remains 
constant by varying U. Besides these peaks, one also 
observes peaks in correspondence to the transition tem- 
peratures separating homogeneous and CO phases. Of 
course, in the case of a reentrant phase, one correspond- 
ingly observes two peaks relative to the transitions. The 
behavior of the specific heat at n = 0.5 as a function of 
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FIG. 13: (a) The specific heat as a function of the temperature 
for z = 3, V = 1, n = 0.5 and different regions of U: (a) 
attractive interaction U = —2, . . . , —0.3; (b) U negative and 
close to zero; (c) U positive and close to zero; (d) U large and 
positive. 



the temperature is shown in Figs. 13a-d. For U nega- 
tive and large (U = —3 V, — 2V), there is no CO and C 
exhibits only one peak at high temperature at position 
Ti w 1; by increasing U (U = -V, -0.5V, -0.3V) a sec- 
ond peak appears at low temperatures T 2 « 0.1 and the 
position of the first peak decreases (see Fig. 13a). Fur- 
ther increasing U, a reentrant phase transition from the 
homogeneous phase to the CO state occurs; correspond- 
ingly, two new peaks are observed around the two critical 
temperatures T 3 w 0.4 and T 4 « 0.15 (see Fig. 13b). For 
attractive on-site interactions the possible excitations are 
creation and annihilation of doubly occupied states, and 
the charge peaks at Ti and T 2 are mainly induced by 
ip^; these peaks move towards low temperatures, with 
T 2 — > as U — > 0. For U small and positive and close to 
zero, there is a phase transition but no reentrant phase, 
T 3 = 0. The specific heat exhibits a three-peak structure 
(see Fig. 13c): one peak at high temperature (Ti w 0.7), 
a second peak at low temperatures (T 2 « 0.02), and a 
third peak around T 3 w 0.4. For U positive and large 
(see Fig. 13d) only one peak at T 3 w 0.5 appears due to 
the phase transition to the CO state. The behavior of 
the specific heat at n = 0.75 as a function of the tem- 
perature is shown in Figs. 14a-c. For n = 0.75 there is 
no reentrant behavior. At low temperatures the system 
is in a CO state and exhibits a phase transition at T c 
to a homogeneous phase. For U negative and large only 
one peak appears at T 3 = T c , due to the phase transi- 
tion. For U negative and small, besides the peak at T 3 , 
one observes another peak T 2 at low temperatures which 
vanishes as U approaches zero. For < U < 3 V, there 
is only one peak at T 3 = T c ; contrarily to the case of 
n = 0.5, the peak T 2 is not observed for U positive. For 
U > 3V (z = 3) there is no phase transition and one 
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FIG. 14: The specific heat as a function of the temperature 
for z = 3, V = 1, n = 0.75 and different regions of U: (a) 
large attractive interaction U = —2,...,—l; (b) U negative 
and close to zero; (c) U large and positive. 



broad peak of much more intensity is observed. The be- 
havior of the specific heat at half filling as a function of 
the temperature is shown in Figs. 15a-b. For n = 1, a 
CO phase is observed for all values of U < Uq. As it is 
shown in Fig. 15a, the specific heat exhibits one peak 
situated at T 3 = T c moving towards low temperatures as 
U increases. At T = 0, for U < Uq, one sublattice (A) 
is empty, whereas the other sublattice (B) has all sites 
doubly occupied; there are neither singly occupied sites 
nor neighbor sites occupied (A^ = 0). The energy of 
the ground state is non-degenerate. For U > Uq, one ob- 
serves a homogeneous state: all sites are singly occupied 
and the energy of the ground state is infinitely degener- 
ate since the spins are arbitrarily oriented. At T = and 
U = Uq, there is a phase transition from a non-degenerate 
level to a degenerate one; therefore one expects a discon- 
tinuity in the internal energy E and, correspondingly, a 
divergence in the specific heat in the limit U — > Uo- This 
is clearly observed in Fig. 15b. 

Other important quantities, useful for studying the 
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FIG. 15: The specific heat as a function of the temperature for 
z — 3, V — 1, n — 1 and different regions of U: (a) values of 
t/ far from the critical point Uo; (b) values of U approaching 
U . 

critical behavior of the system, are the charge and spin 
susceptibilities. In fact, anomalies in their behaviors 
clearly signal the onset of a CO state. The charge suscep- 
tibility can be defined separately for the two sublattices: 

A _ dn A B _ dn B 

Xc ~ cV ' Xc ~ 8ft ' 

where the total charge susceptibility is given by Xc = 
(Xc + Xc ) /2- By studying separately Xc an d Xc > & is 
manifest their different contribution to the total charge 
susceptibility. By varying the particle density at fixed 
(low) temperature, one finds a critical value of the parti- 
cle density (n c ) at which the susceptibilities Xc an d Xc 
show a discontinuity: below n c , one finds = Xc = Xc, 
whereas above n c one observes the depletion of sublatticc 
A and the corresponding accumulation of the particles in 
the sublattice B. For n > n c , the susceptibility xt be- 
comes negative, since n A is a decreasing function of n. 
As a result, the system is in a CO phase with "almost 
empty" shells separated by filled shells. Of course, at 
half filling and in the limit T — ► 0, one finds ha = and 
tib = 2 and the susceptibility vanishes. In Fig. 16, we 
plot the charge susceptibility as a function of the particle 
density for a given value of U and several values of the 
temperature. At low temperatures, \c shows a quasi-one 
or -two lobe structure, depending on the value of the on- 
site potential. For z = 2, i.e., a ID chain, the lobes have 
a regular shape [1]. For z > 3, the curvature of the lobes 
is different above and below n c and the observed disconti- 



FIG. 16: The charge susceptibility as a function of the particle 
density for z = 3, V = 1, different values of the temperatures 
and (a) U = -1; (b) U = 2.1. 

nuity signals the transition to a phase where translational 
invariance is broken. The vanishing of \c for T — > at 
quarter and half filling indicates that a full CO state is es- 
tablished. In Figs. 17a-b, we plot Xc as a function of the 
temperature, for U = V and U = 2.1V, and for different 
values of the filling (n = 0.25, 0.5, 0.75, 1). For n = 0.25 
there is no charge ordering and one finds x^ = Xc = Xc- 
At n = 0.5 and n — 0.75, the system is in a CO phase 
characterized by a different filling of the shells for T < T c : 
as a consequence, xt Xc ■> an d Xc presents a disconti- 
nuity. At half filling, there is charge ordering but the 
susceptibility is the same in the two sublattices, also for 
T < T c . The limit T — > of the charge susceptibility dra- 
matically depends on the value of the particle density: Xc 
is finite for values of the particle density corresponding to 
translational invariant states, whereas it decreases with 
T and vanishes for fillings corresponding to CO states 
characterized by alternating empty and filled shells. In 
the limit of high temperatures, the charge susceptibility 
tends to a constant value which does not depend on U 
and z but only on n according to the law [f] 

lim Xc = a(n), (55) 

T — >oo 

where a(n) = n(2 — n)/2. 

The spin magnetic susceptibility Xs in zero field, given 
in Eq. (41), can also be defined separately for the two 
sublattices. In Figs. 18a-b we plot the spin susceptibil- 
ities Xs, xf an d xf as functions of the temperature for 
two representative values of U (U = —V and U = 2.1V, 
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FIG. 17: The charge susceptibility \ c as a function of the 
temperature for z = 3, V = 1, n = 0.25, 0.5, 0.75, 1 and (a) 
[/ = 1, (b) U = 2.1. 



FIG. 18: (a) The spin susceptibilities Xsi Xs 1 an d X? 85 a 
function of the temperature for V = 1, n = 0.25, 0.5, 0.75, 1 
and for (a) [/ = -1 and (b) U = 2.1. 



respectively) and for different values of the filling. For 
U = —V, the spin susceptibility vanishes at zero tem- 
perature for all values of the filling: all the electrons are 
paired and no alignment of the spin is possible. By in- 
creasing T, the thermal excitations break some of the 
doublons and a small magnetic field may induce a fi- 
nite magnetization: Xs augments by increasing T up to 
a maximum, which might be different for the two sub- 
lattices, then decreases. For n = 0.25 and n — 0.5 the 
system is not charge ordered and Xs has the same value in 
the two sublattices. Conversely, for n = 0.75 and n = 0.9 
the system is charge ordered when T < T c : Xs assumes 
different values in the two sublattices. 

For U = 2.1V, when the system is not charge ordered, 
xf = xf — Xs diverges for T — > 0, whereas for CO states 
xf vanishes and xf diverges in the same limit. At low 
temperatures, only the electrons belonging to sublatticc 
A are paired and, as a consequence, xf — 0- At half 
filling, for both U — —V and U — 2.1V, even in the 
presence of charge ordering, Xs has the same value in 
the two sublattices A and B. It is easy to check that, 
for high temperatures, the spin susceptibility decreases 
with the Curie law: limT^oo Xs = a(n)/T in the entire 
(U, n) plane, where a(n) is the same z and U independent 
function appearing in Eq. (55) [1]. As a consequence, 
in the limit of high temperatures, the ratio Xc/Xs is an 
universal function of T, namely: lim^— oo (Xc/Xs) = T. 

To conclude this section, we report some results ob- 
tained for the entropy as a function of T, and U, showing 



that also this quantity is a good indicator of the onset of 
a CO state. The standard way to compute the entropy 
is by integrating via the integral of the specific heat: 

S(T) = S(0) + J Q T ^dT. 

However, as extensively commented on in Rcf. [1], this 
expression can not be easily handled since it requires the 
calculation of S(0), which is generally not an easy task. 
A more convenient formula is given by: 

In our framework of calculations, we can readily deal with 
this expression since it requires only the knowledge of the 
chemical potential. In Figs. 19a and 19b we plot the en- 
tropy as a function of the temperature for n = 0.5 and 
n = 1, respectively. One can notice an abrupt change 
of the entropy curves when the critical temperature is 
reached. One may also notice that, in the limit T — > 0, 
S(0) is finite for n — 0.5 due to the degeneracy of the 
ground state energy. On the other hand, 5(0) is zero at 
half filling because there is no degeneracy of the ground 
state. The U dependence of the entropy is rather dra- 
matic in the neighborhood of the values at which a zero 
temperature transition occurs, as it is evident from Fig. 
20. At low temperatures, the entropy presents a step-like 
behavior and becomes rather insensitive to variations in 
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FIG. 19: The entropy S as a function of the temperature for 
2 = 3, V = 1, U = -1, 1, 2.1 and (a) n = 0.5, (b) n = 1. 

U for sufficiently large on-site repulsive and attractive 
interactions. For < n < 1/z, one observes an increase 
of the entropy at U = in the limit T — > 0: even if 
there is no transition, there is a change of the distribu- 
tion of the electrons from a configuration characterized 
by doubly occupied sites to one with singly occupied sites, 
which is less ordered. When 1/z < n < 2/z, the entropy 
presents a discontinuity around U — (due to the phase 
transition) which becomes more pronounced as the tem- 
perature decreases. In the region 2/z < n < 1, at low 
temperatures, one finds a rather sharp increase of the 
entropy near Uo(n), where the system undergoes a tran- 
sition from a CO phase to a translational invariant state, 
less ordered. 



V. CONCLUDING REMARKS 

Statistical models on the Bcthe lattice are of consid- 
erable interest since they admit a direct analytical ap- 
proach for a number of problems that may be otherwise 
intractable on Euclidean lattices. In this paper, we have 
evidenced how the use of the Green's function and equa- 
tions of motion formalism leads to the exact solution of 
the extended Hubbard model on the Bethe lattice in the 
narrow-band limit. We provided a comprehensive and 
systematic analysis of the model by considering relevant 
thermodynamic quantities in the whole space of the pa- 
rameters n, T, U and V and we obtained the finite tem- 
perature phase diagram, for both attractive and repulsive 



(a) (b) 




FIG. 20: The entropy S as a function of U for z = 3, V = 1, 
for temperatures varying in the interval (0,1) and (a): n = 
0.25, (b): n = 0.5, (c): n = 0.75, (d): n = 1. 



on-site and intersite interactions. 

The phase diagram dramatically depends on the sign of 
the nearest-neighbor interaction. In the attractive case, 
there is a critical temperature T c , located by the diver- 
gence of the charge susceptibility, at which there is a 
transition from a thermodynamically stable to an unsta- 
ble phase; the latter is characterized by phase separation. 
The surface separating the two phases has a rounded 
vault in the 3D (n, U, T) space, which splits in two for 
U > U c {z). Below T c , the same critical value of the on- 
site potential separates the two possible configurations 
of phase separated states, namely clusters of singly or 
doubly occupied sites. For repulsive nearest-neighbor in- 
teractions, the phase diagram has a richer structure: we 
found a transition temperature below which translational 
invariance is broken. The Bethe lattice effectively splits 
in two sublattices with different thermodynamic proper- 
ties. As a result, a charge ordered phase, characterized 
by a different distribution of the electrons in alternat- 
ing shells, is established for n > 1/z. The CO phase 
is energetically favored as demonstrated by the study of 
the Helmholtz free energy. By investigating the behavior 
of several thermodynamic quantities, we attained a full 
comprehension of the phase diagram. The study of the 
particle density and of the double occupancy is particu- 
larly enlightening to unveil the distribution of the parti- 
cles on the sites of the Bethe lattice. The specific heat 
exhibits not only high and low temperature peaks due 
to charge excitations induced by the Hubbard operators 
and tj}^\ but also peaks due to the phase transi- 
tion from the homogeneous phase to the CO state (or 
viceversa). By studying separately the contribution to 
the charge and spin susceptibilities coming from the two 
sublattices, the onset of a CO phase is also signalled by 
the separation of the values of the sublattices quantities 
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X^ s and Xcs ■ The different population of the two sublat- 
tices in the CO phase is reflected by the divergence or the 
vanishing, in the limit T — ► 0, of the spin susceptibilities 
xf and Xs- If the electrons are paired, no alignment of 
the spin is possible (xf — ► 0), whereas single occupation 
of the sites leads to an alignment (xf — ► oo) even when 
the magnetic field is turned off. 
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APPENDIX A: The extended Hubbard model on 
the Bethe lattice in the presence of an external 
magnetic field 

In the presence of an external magnetic field h, the 
Hamiltonian of the AL-EHM reads: 



H =-^Y,n{i) + uY,D{i) 

i i 



(Al) 



where n 3 (i) is the third component of the spin density 
operator 

n 3 (i) = n T (i) - n^i) = c{(i)c T (i) - c\(i)ci(i). (A2) 



The formulation given in Section II for the case of a Bethe 
lattice must be generalized in order to take into account 
the breaking of rotational invariance in the spin space: 
one has to distinguish the two components - in the spino- 
rial notation - of the fermionic fields. By exploiting the 
decomposition (25), it is not difficult to show that the 
presence of the magnetic field affects only Hq. In this 
representation the equations of motion of the Hubbard 
operators become 



[&(i),H ] =-(fi + ah)^(i), 
[r] a (i),H Q ] =-((!, + ah - U) r] a (i) 

and Eqs. (38) are modified as: 

e 0(n+2h) , e 0{2n+h-U) 



(A3) 



(ni(*)>o = 
(D(i))o = 



e f3h + e (3n + e /3(n+2h) + e /3(2n+h-U) ' 

e /3 M + e f3(2n+h-U) 
e f3h _|_ e (3n + e 0(p.+2h) + e /3(2p,+h-U) ' 



o 0(2n+h-U) 



(A4) 



e 0h _|_ e 0fj, _|_ e 0(n+2h) _|_ e 0(2^+h-U) ■ 
e^(e 2 ^ h - 1) 

e /3h ZjZ e 0n _|_ e /3(/i+2?i) _|_ e 0(2n+h-U) ' 



All the rest of formulation developed in Sections II, III 
and IV follows easily. It is worth noticing that the mag- 
netization takes the simple expression: 

(n 3 (*)) = tanh(/3/ l )[(n(i)) - 2{D(i))]. (A5) 
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